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ABSTRACT 

We investigate the use of global demons, a 'canonical dynamics', as an 
approach to simulating lattice regularized field theories. This determinis- 
tically chaotic dynamics is non-local and non-Hamiltonian, and preserves 
the canonical measure rather than 5(H — E). We apply this inexact dy- 
namics to the 2D XY model, comparing to various implementations of 
hybrid Monte Carlo, focusing on critical exponents and critical slowing 
down. In addition, we discuss a scheme for making energy non- conserving 
dynamical algorithms exact without the use of a Metropolis hit. 
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1. Introduction 



One traditional dynamical approach to simulating ensemble averages has been molec- 
ular dynamics (MD) algorithms. In the simplest of these, micro-canonical simulations, 
conjugate momenta are introduced for each degree of freedom in the ensemble, and the 
resulting system is time-evolved according to Hamilton's equations of motion. The re- 
versibility of Hamiltonian evolution then ensures detailed balance, i.e. that the simulation 
is a Markov process. If the system is sufficiently large, and the interactions sufficiently 
complex, one usually imposes the quasi-ergodic hypothesis, and hopes that the simulation 
will explore the desired ensemble. Unfortunately, it is known that Hamiltonian evolution 
conserves energy, and is therefore not ergodic. In fact, microcanonical simulations intro- 
duce an explicit factor of 8{H — E) into the measure of the ensemble being simulated. In 
order to use an MD algorithm to obtain the correct ensemble, some additional method 
must be introduced to integrate over the different energy surfaces E. 

One method of dealing with this difficulty is embodied in the hybrid molecular dynam- 
ics (HMD) and hybrid Monte Carlo (HMC) algorithms [pJ-0. In these algorithms, the 
micro-canonical equations of motion are integrated along a 'trajectory' for a time T, after 
which the momenta are touched with a heat bath, changing the energy of the system. As 
with any numerical integration of Hamilton's equations, finite step size errors will build 
up along the micro-canonical trajectories, leading to systematic errors in the ensemble 
generated by HMD. Although these errors can be controlled by making the step size suffi- 
ciently small, this can become costly. The HMC algorithm is designed to correct for these 
dt errors. It does this by treating the configuration at the end of a micro-canonical tra- 
jectory as a proposal for a global update of the system, which is then accepted or rejected 
according to a Metropolis hit. If the equations of motion are reversible, this sequence of 
configurations is then a Markov chain which, given ergodicity, is guaranteed to produce 
the correct ensemble. The HMC and HMD algorithms are currently widely used in lattice 
gauge theory simulations, especially in systems involving dynamical fermions. Typically, 
HMC is used in theories where the action can be expressed as the volume integral of a 
local function, while HMD is used when it cannot. 

Although these algorithms are generally quite robust, they do have one weakness, 
related to critical slowing down. Associated with any observable O is an autocorrelation 
'time' To, the time scale required for the simulation to produce a statistically independent 
measurement of O. This autocorrelation time will generally depend on the correlation 
length of the system as a power law[Q, 

T = AC. (1) 

where z is the dynamical critical exponent, and A and z will depend on O. Critical slow- 
ing down occurs whenever z > 0. To represents the typical amount of simulation time 
it takes for a local change in O to propagate across a correlated cluster. For multi-scale 
algorithms, such as cluster algorithms, one can hope to obtain z = 0, since these algo- 
rithms are designed to change an entire correlated cluster simultaneously. Unfortunately, 
these algorithms are not easily generalized from one model to another, and have not yet 
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been implemented for most lattice gauge theory models. At the other extreme, it is dan- 
gerously easy to obtain z = 2 with an algorithm which involves only local updates. This 
just corresponds to diffusive transport of fluctuations through the correlated cluster. In 
principle, one should be able to do no better than z — 1 with a local algorithm, since the 
local algorithm will restrict fluctuations to a finite propagation speed^j. One can obtain 
z = 1 [f|,|6[] in HMC and HMD, but to do so requires a correlation length dependant tuning 
of the trajectory length. If the trajectories are too short, then the frequent randomiza- 
tion of the momenta causes the motion of the system through phase space to be diffusive 
(resulting in a z of 2), while if the trajectories are too long the quality of statistics goes 
down because energy conservation correlates the measurements along each trajectory. In 
addition, the optimum trajectory length is likely to differ for different observables, forcing 
an inefficient trajectory length for some of them. (This problem is especially severe for 
HMC, where only one measurement is allowed per trajectory. Running with trajectories 
twice as long as is necessary is therefore equivalent to a factor of two slow-down in the 
code.) 

We would like to contrast these hybrid algorithms based on micro-canonical evolution 
to a purely dynamical approach we term global demons. Like the hybrid algorithms, 
global demons are an easily implemented, very general approach to simulating ensemble 
averages. Unlike the hybrid algorithms, the global demon formulation has nothing to do 
with Hamiltonian dynamics. For example, it can be defined completely in terms of coor- 
dinates alone if so desired: the presence or absence of a symplectic structure is irrelevant. 
This should be contrasted to the MD based algorithms, where a partition function whose 
action depends only on coordinates is usually augmented to include fictitious conjugate 
momenta in order to define a Hamiltonian or Poisson structure. Another way to say this 
is that while evolution using Hamiltonian dynamics generates a microcanonical ensemble, 
evolution in the global demon system generates a canonical ensemble. The global demon 
equations of motion are deterministic and time reversal invariant, and are designed to 
evolve through the physically accessible regions of configuration space (it is not a phase 
space) such that the trajectory fills configuration space with a density reproducing the 
correct ensemble. Consider, as an illustration, 1000 points in a phase space (q,p) which 
lie equally spaced on a unit circle, as shown in the left column of Fig. 1, at t — 0. The 
points are connected in order to see how neighboring points behave. Under microcanoni- 
cal evolution of the 1-d harmonic oscillator equations of motion (t = 1 is the natural time 
scale of the dynamics), 

q = p, p = -q, (2) 

this circle will be preserved, and the points will rotate, preserving the figure at all later 
times. In contrast to Hamilton's dynamics, global demon dynamics will result in a rapid 
dissemination of neighboring points through the space. The time evolution of the circle for 
a 1-d harmonic oscillator Hamiltonian H = (p 2 + q 2 )/2 is shown at time t = 10 and t = 20, 
in our canonical dynamics. What is striking is the speed at which neighboring points on 
the circle evolve to opposite sides of the the space. Even at t — 20, one can see that the 
phase space density is nearing the desired ensemble exp(—j3H). In the right column of 
Fig. 1, we have the same situation for an SU(2) Hamiltonian H = J 2 /2, whose phase 
space, parameterized by (J x , J y , J z ), is the unit sphere. Here an initial condition of a circle 
at J z = 0.5 is evolved in a similar manner. Again, the rapid divergence of neighboring 
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points is striking. Since global demons generate a deterministic chaotic dynamics, the 
danger of diffusive motion through phase space present in the hybrid and other stochastic 
algorithms is absent here. The one drawback to the global demons approach is that we 
have been unable to determine a way of making it exact, i.e. removing the dt errors in 
the ensemble arising from finite step size. 

In this article, we investigate an application of the global demon algorithm to a lattice 
field theory. We are interested in understanding the critical properties of the dynamics 
near phase transitions, and how tuning the dynamics can improve convergence. We also 
want to examine the correctness of measurements in this inexact dynamics, as well as ways 
to make it exact without recourse to stochastic techniques. We choose the 2D XY model, 
since it has been well studied in the past in a variety of algorithms, and has an infinite 
order phase transition around which we can investigate the manifestations of critical 
slowing down. In section 2, we present the global demon dynamics for unconstrained 
systems. The implementation of this dynamics for the XY model in section 3 explores the 
behavior of simulations as various parameters of the algorithm are tuned. We find that 
the algorithm is quite robust, obtaining good results with little tuning. In addition, we 
compare the critical behavior of the algorithm to several implementations of HMC We 
conclude in section 4. Finally, in an appendix, we present a scheme which, in principle, 
should remove the finite dt errors from the algorithm dynamically. Although our present 
implementation of the global demon algorithm is not exact, we have chosen to compare 
our results to HMC, rather than HMD. The main reason for this choice is that we are 
treating the HMC results as a control, and would like them to be as free of systematic 
errors as possible. Since the HMC and HMD algorithms are so similar, however, it is 
quite likely that qualitative conclusions about the critical behavior of HMC will also be 
correct for HMD. 



2. Global Demon Dynamics 



Let us consider a system characterized by an action S(x) and coordinates x — (x±, ...,x n ). 
The ensemble averages of this system will have the generic form 

(O) = \ J Vn(x) e-WO, (3) 

where 

Z = f Vn{x)e~^ x \ (4) 

is the normalization, and the measure Vfi(x) might include constraints (for example, 
symmetries associated with a Lie algebra). In this article, we are only concerned with the 
situation when the measure is trivial, T>x. If the variables x = (q,p) include canonically 
conjugate coordinates and momenta, S(x) can be taken as a Hamiltonian: H(p, q) = S(x). 
Otherwise, in what is more or less standard practice, conjugate momenta are introduced 
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and added to S(x) to make the exponent in Eq. (|3]) resemble a Hamiltonian 



1 n 

H(p,x) = -J2Pi+S(x), (5) 

i=i 

and the measure is modified to VxT>pexp(—f3H), with appropriate normalization. Molec- 
ular dynamics is now easily implemented, leading to the equation of motion 

Xi = Pi, Pi = 5 , (i=l,...,n). (6) 

OXi 

This dynamics, combined with momentum refreshes and global metropolis hits, produce 
the HMD and HMC algorithms. 

Let us now pass to the canonical dynamics for such a system |7|-||. In contrast 
to the microcanonical dynamics, the energy and the symplectic structure are no longer 
preserved. Rather, the measure itself is preserved directly by the dynamics. We can define 
such dynamics in many ways. For instance, instead of Eq. @, we could take 

(3 dG(w) 

Xi = -K - Fi(x) , (i = l,...,n). (7a) 

n dw 

Here G(w) and Fi(x) are arbitrary functions of a global demon variable w and coordinates 
x, respectively, and k is a coupling constant. The number of global demon interactions (the 
right hand side of (7a)) is unrestricted. In this example we have used 1, while 2 or 3 are 
usually sufficient, regardless of n. This type of treatment can be viewed as a deterministic 
version of Parisi and Wu's stochastic quantization [0, ||. An alternative formulation of 
the x dynamics which includes a relic of the underlying Hamiltonian H(p,x) is 

ii = Pi- — F u (x), (7b) 
n aw i 

dS(x) k 2 /3 dG 2 „ , v v 

Vi = 7^ -I— F 2i {p). (7c) 

OXi n dw 2 

An important observation here is that while we can retain a Hamiltonian sub-structure to 
the dynamics, it is not responsible for the ergodicity in the full configuration space, and can 
be retained or altogether removed. This will have some effect on the convergence, since the 
Hamiltonian forces can provide additional decorrelation. In non-equilibrium simulations, 
it is more convenient to use (7b-c) since there is a closer link to the thermodynamics of 
H(p, x)PT|. Eqs. (7b-c) also have a microcanonical limit when K a = 0. Since we are 
going to compare the global demon approach to HMC, we retain the momenta to have a 
greater parity between the two algorithms. 

With the introduction of the global demons w a , a larger configuration space {(f)} must 
be defined, where <p = (xi, ...,p n , u>i, .■■w m ) . In 0— space we can define a new action /, 
which is determined by the equations of motion (7) in the following way: 

-in m 

f(x,p,w) = S(x) + -Y.pl + E G «M, P f = e-M. (8) 

1 i=\ a=l 
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Unlike H(p, x) in molecular dynamics, / is not preserved by the equations of motion. 
While the definition of / is not unique (in the sense that the measure for the variables w 
is arbitrary), it is natural as well as convenient in determining the dynamics of the global 
demons, by providing a particular solution to the continuity equation below. Eq. (H) de- 
fines the 0— space measure as DxDpDw exp(— /?/) = pf T><p. The global demon dynamics 
can then be determined by requiring that pj be a stationary solution of a generalized 
Liouville (continuity) equation in configuration space: 

=^ (9 ) 

at d< Pi 

This is equivalent to requiring that the master equation, enforcing conservation of prob- 
ability under evolution of the ensemble, be satisfied. 

The equations of motion for the demons are now found by requiring that they, com- 
bined with the generalized dynamics of (7), satisfy Eq. @. A direct substitution of pf 
and <pi (Eqs. (7b-c) ) into Eq. (|9|) allows one to solve for w: 

«! / OS dF lt \ 
«>i = — W^-^T h ( 10 ) 



«1 


[t>F™ 

\ OXi 


dF lt 


n 


dxi 


n 


\0F 2iPi - 


dF 2l \ 
dpi J 



w 2 -- 

If we had chosen to neglect the momenta and used (7a), the form of Eq. (|T0|) would be 
unchanged. Eqs. (7) and ( |T0D define a dynamics which by construction preserves the 
measure Eq. @. (It is worth noting that while we have taken an exponential form for the 
density, in general we can take an arbitrary function p and still use this same procedure.) 

Microcanonical dynamics preserve the phase space volume exactly, since the divergence 
of the equations of motion, 

d A + d A (in 

dqi dpi ' 

trivially vanishes by Hamilton's equations of motion. The global demon equations of 
motion (7), (fL0|), on the other hand, do allow for fluctuations in the 0-space volume, which 
can be quite large. Writing these equations as (pi = J-i{(p), the divergence is explicitly 

dip, (3 ( dG 1 dF li dG 2 dF 2l \ 

wr-n{ Ki d^^x- +K ^ 2 ^p-)- (12) 

This local 'breathing' of <p— space is controlled by the arbitrary functions G and F. Al- 
though this behavior is not microcanonical, there is nevertheless an invariant quantity, 
called the pseudoenergy £, which is preserved: 

e = f(x, P ,w) + \ fdt 1 ^. (is) 



(3 Jo 

One can check directly that £ = 0. 
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There is clearly some freedom in defining the dynamics: the functions G and Fi and 
the coupling strength k. The only restriction on G(w) is that the measure Eq. (|8]) is nor- 
malizable; in general the auxiliary variables w can have any desired measure. In practice, 
highly non-linear functions are impractical since they will require small integration time 
steps. For these reasons, it is convenient to take G = w 2 /2 or G = w A /4. A necessary 
condition for Fi is for it to be at least linear in its argument, the minimal requirement for 
the existence of the fluctuations in the volume (|T2"D . The precise relation to the fluctua- 
tions in a volume V, or equivalently, the instantaneous 0— space compressibility, can be 
found using the divergence theorem 

dV /3 r ( dG 1 dF lt dG 2 dF 2t \ 

In this paper, we do not explore the effect of different choices of Gi, F^ Such studies have 
been done on smaller systems [pi -H. 



Finally, we observe that the equations of motion (7), fllOD will have no stable fixed 



pointsfl^j. This is the case since the sum of the Lyapunov exponents is related to the 



average rate of change of total volume of 0— space. By the Liouville equation, this will 
necessarily vanish [pf] : 

n+m Qi 

EA, = <S> = 0. (15) 



3. Implementation for the 2D XY Model 



The 2D XY model consists of spins located on the sites of a two dimensional square 
lattice, which are free to rotate in the plane. The action is given by 

v{6) = ~J2 ReU i u j = - E cos & - j)> ( 16 ) 

<ij> <ij> 

where the sum is over nearest neighbors, and the Uj = e are elements of U(l) located 
at each lattice site j. In two dimensions, this model exhibits a Kosterliz-Thouless phase 



transition near (3 ~ 1 [ 14 —[V? . Above the phase transition, the dynamics is dominated by 
dissociated vortex-antivortex pairs. These pairs become tightly bound below the phase 
transition, where the dynamics is dominated by spin waves. The K-T phase transition is 
infinite order, characterized by an exponentially diverging correlation length (£): 

^a^MkiT-nm- (17) 

Numerical simulations indicate similar critical behavior for finite lattices^, O. Near 
the critical temperature, the system experiences an exponential increase in the correlation 
length £, which can lead to critical slowing down in simulations by virtue of Eq. (p. 
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3.1 Equations of Motion 



The implementation of global demons to the XY model is straight forward. Following 
the equations of motion (7b-c) and fllOD , we have 



0* 



Pi 



Pi 



K 2/3 .3 

-W2 sm Ui 



n 

dV{9) 



l,...,n) 



Ktp 3 

-w x pi. 



n 



Fi{6) = sin 3 ^ is chosen to respect the periodicity in 9, while Gi(p) = Pi has no such 
restriction. This choice was motivated only by simplicity, and in general, we could take 
more complicated interactions, and include additional global demons. The corresponding 
equations for the global demons are then 



Wi 



w 2 



n 



n , 



^E^^sin 3 ^-3^sin 3 ^cos^ 



(19) 



We used leapfrog integration, which included a Taylor expansion so that the 0(dt 2 ) errors 
in a time step cancel. The pairs q, W\ and p, W2 were updated in alternate steps. w\ was 
taken with q since they both involve momenta, and W2 with p since they both involve co- 
ordinates. Our general studies of the systematics of the model under tuning of parameters 
were performed on a 16 2 lattice, while the studies of the critical exponents were done on a 
64 2 lattice, to allow longer correlation lengths. The particular choice of functions F leads 
to a small non-ergodicity for this particular system: the momentum zero mode cannot 
change sign. (We have corrected for this by occasionally (every 64 trajectories) refreshing 
the momenta, using the same procedure as in HMD. In general, this is probably a good 
idea, to ensure that the evolution is ergodic. This particular non-ergodicity could also 
have been corrected by a small modification of the equations of motion.) We have veri- 
fied that the equations of motion are correct to 0(dt 3 ) on a time step, leading to 0(dt 2 ) 
systematic errors in observables, by computing the dt behavior of several observables in 
Fig. 2, demonstrating the quadratic behavior of the systematic error. 



3.2 Hybrid Monte Carlo 



We have used the critical properties of HMC as a benchmark for comparison of our 
global demon approach, studying three of its variations^]. The equations of motion are the 
same as those used for global demons, except with Ki — 0. By modifying the length of the 
HMC trajectory between momentum refreshes, we modify the decorrelation time@. The 
first variation, denoted HMC-1, has trajectories of length 1, where the highest frequency 
of the free theory is (2ir)~ 1 . While this 'standard' choice is easy to implement, it suffers 
from severe critical slowing down, with z = 2 in Eq. (|I]). The critical behavior should be 
improved by choosing the trajectory length proportional to the spatial correlation length 
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£ of the system^]. The two variations we consider are denoted HMC-S, for T — £, and 
HMC-L, for T = 27r£. Again we point out that, in order to make this choice, we require 
the very information which we are attempting to measure. We were fortunate to have 
previous results for £ available to usjjrTJ, but in general this is not likely to be the case. 

The integration time step was kept fixed at dt = .1 along trajectories of length T. 
This value of dt was chosen so that the acceptance rate of the global Metropolis hit in 
the HMC algorithm was approximately 80%. For HMC, it is necessary to use random 
trajectory lengths for optimum relaxation ||: we chose T uniformly distributed on the 
interval (.5(T), 1.5(T)], with (T) = 1 for HMC-1, and (T) « f , 2tt£ for HMC-S,L. We also 
made several runs of HMC-L using exponentially distributed random trajectory lengths. 
We found that this does not lead to any improvement over the runs with uniformly 
distributed trajectories, and may even have resulted in slightly noisier measurements. 
One interesting difference between HMC and global demons is that, for global demons, 
T simply denotes time between measurements along a single trajectory - the evolution 
of the simulation is completely unaffected by the choice of T. At the end of a HMC 
trajectory we performed a global metropolis hit, after which we performed measurements 
and refreshed the momenta by choosing new, gaussian distributed pj. 

3.3 Coupling Strength Dependence 

In the micro-canonical algorithms HMC and HMD, the rate at which the simulation 
covers phase space in the non-ergodic directions (i.e. changes energy) is controlled by the 
time between momentum refreshes. If the trajectories are too long, the system changes 
total energy very slowly, leading to autocorrelations on timescales proportional to the 
trajectory length. If the trajectories are too short, on the other hand, the motion between 
energy shells is rapid, but motion in the micro-canonical directions is diffusive, leading 
to a dynamical critical exponent of 2. This means that, for large correlation lengths, 
the efficiency of the algorithm can vary as a power of T/T opt , where T opt is the optimum 
trajectory length. As shown in Refs. f|], £§], T opt should be proportional to the correlation 
length of the system, which is not known a priori. Thus it is necessary to perform a 
sensitive tuning which depends on a parameter measured in the simulation. 

In the global demon algorithm, on the other hand, the parameter which controls energy 
(or action) non-conservation is just the coupling k of the demons to the system. In the 
limit k — > 0, the demons decouple and ergodicity is lost. If k becomes too large, the 
equations of motion will suffer the characteristic instabilities of discretized dynamics. In 
contrast to the HMC and HMD algorithms, however, we can make an a priori choice of 
k which works quite well at all values of the correlation length. To do this, consider the 
change in total action, AS = S(<p2) — S(<j>i), in a single time step, 

(99 

AS ~ At -^—<pi , (20) 

d(f)i 

where S now refers to the total global demon action / in Eq. (§). Using equations (7)-([10|) 
and (POD, we see that < \dS/dt\ > is proportional to k, with a constant of proportionality 
of order one. Thus, the change in S along a trajetory of length T should be AS ~ kT. To 
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set the scale of AS, we can compute its expectation value if two consecutive measurements 
are totally decorrelated: 



a = yf(A!P) = y/2((S(<f>y) - (W) = i , (21) 

where C s is the specific heat of the system. In order for the action to decorrelate between 
measurements (with T = 1), we conclude that the optimal choice of k is likely to be 
k ~ O(yfn). Note that this philosophy of forcing large fluctuations in energy along a 
trajectory is inherently different from HMC, where (in order to avoid prohibitively low 
acceptance rates) AS" is 0(1). 

To investigate the behavior of the global demons algorithm under tuning, we ran a 
series of simulations on a 16 2 lattice. (Except where otherwise indicated, we used an 
integration time step of dt = 0.1 and measuring at intervals T = 1 = lOdt; we will call 
this time between measurements a trajectory length, even though no momentum refresh 
is performed). Along the trajectory, we determine the square of the change in action 
between measurements, AS 2 = [S(t) — S(t — T)] 2 , which is then averaged along the entire 



trajectory to obtain y(AS 2 ). The result gives a guide as to how fast the trajectory can 
diffuse through configuration space. In Fig. 3, we plot the quantity AS, which we call the 
diffusiveness, defined by 

(AS 2 ) 



AS = ^- -, (22) 

a 

as a function of coupling strength k, for simulations at representative values of (3 both 
above and below the phase transition. In the limit = 0, the equations of motion 
(p!8|)-(]T9|) are microcanonical and S is preserved, as indicated in the figure. For small 
couplings, the microcanonical component of the dynamics is only slightly perturbed by 
the canonical component, and the ergodicity is weak. The convergence times here are 



quite large pq|. For Ki ~ y2n ~ 23, the value of AS can be seen to saturate near unity, 
the value expected when two consecutive measurements are uncorrelated; here the steps 
are quite large through configuration space. For larger values of the coupling, the change 
in action remains saturated. Convergence is also generally slower for larger k, since the 
additional decorrelation produced by the microcanonical component to the dynamics is 
reduced. The reduction of AS in Fig. 3 for (3 ~ 1 can be attributed to critical slowing 
down. However, while the correlation lengths become quite large, the dip in AS" is not so 
noticeable. In this respect, critical slowing down does not seem to hinder the dynamics, 
nor does it require any special tuning of k. 

In Fig. 4, the (3 dependence of AS is plotted for simulations at fixed coupling strengths. 
By selecting the k = 32 curve, for instance, we see that we can study both the low and 
high temperature properties of the XY model, as well as the phase transition, without 
modifying k. While there is a small dip in the curves near /3 ~ 1, critical slowing down 
does not seem to strongly effect this measure of the dynamics as one approaches the phase 
transition from either side. An important result is that the couplings k are essentially 
independent of (3 as well as the details of the physics of the model under study. 

It should be emphasized that the runs with k <C \^2n converge slower as k decreases, 
and ultimately do not converge for the microcanonical limit k = 0. We have checked 
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the convergence of the dynamics to the proper ensemble by measuring and subsequently 
histogramming w a and pi along the trajectory and comparing them to their exact analytic 
distributions, finding that convergence is best above k ~ \f2n. We can also examine the 
k dependence of measurements at fixed /3, illustrated in Fig. 5 for (3 = 1/1.1. What we 
see is that the measurements for coupling strengths roughly 5-8 times saturation do not 
exhibit any systematic deviation as k increases. For much larger k (at fixed dt) , there will 
be the characteristic instabilities associated with difference equations. However, the value 
K = y/2n clearly is not near this instability limit, and we can safely use it. In the low k 
limit, we are close to micro-canonical dynamics, and S begins to become approximately 
conserved. This slow diffusion results in long time correlations and poor statistics. This 
figure is typical of other temperatures, above and below the phase transition. 

3.4 Choice of Trajectory Length (i.e. Measurement Frequency) 

An indication of how rapidly measurements decorrelate is shown in Fig. 6. There we 
plot the diffusiveness AS as a function of trajectory length T, for K\ — «2 = 16 and (3 = 1. 
We observe a saturation in the trajectory length near T = 2. Because measurements do 
not effect the time evolution of the global demon trajectory (they are not associated with 
any Metropolis hit or momentum refresh), the choice of frequency of measurements is 
governed by the relative costs of the time evolution and measurement routines. Because 
our measurement algorithm was relatively inexpensive, we chose to use T = 1. 



3.5 Observables 



The observables we measured include the energy E, lattice magnetization M, topo- 
logical charge Q, defined by 

E = ReU t Uj, 

n <ir> 

M = (23) 
n 

Q = -H<? P - 

The sum in Q indicates the sum over all plaquettes of the number of positive topological 
charges occupying that plaquette. (For an exact definition of the topological charge, 
see e.g.Rl9|.) Corresponding to these observables, we can define the specific heat and 
susceptibilities: 

C v = (3 2 n((E 2 )-(E) 2 ), 

X Q = n{(Q 2 ) - (Q) 2 ), (24) 
X m = n(((ReM) 2 ) + ((ImM) 2 )). 

In both our global demon and HMC runs, we started with about 20000/T trajectories 
for thermalization, followed by 160,000/T trajectories of data, where T is the trajectory 
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length. Statistical errors in observables were obtained by binning measurements in bins 
of size 2 n . The errors quoted in our tables use the smallest bin larger than 8tm, where tm 
is the integrated autocorrelation time of the total magnetization. The errors in the sus- 
ceptibilities were obtained from the errors in the corresponding observables by assuming 
gaussian fluctuations on a timescale r , where O is the appropriate observable. A selec- 
tion of our observables are indicated in Table 1. We find that the global demon results 
usually agree with the HMC results within a few a, which indicates that the systematic 
errors are not large. They could, of course, be further reduced by extrapolating to dt = 0. 

3.6 Autocorrelation Functions and Decorrelation Times 

Because we have a dynamical algorithm, the trajectory has a memory, which will 
be reflected in the auto-correlation functions. This can be analyzed by examining auto- 
correlation functions of the observables M, E, Q, and S. When the couplings k are small, 
the dynamics is near the microcanonical limit, and decorrelation is very poor. Typical 
auto-correlation functions for Ki = k 2 = 1 are shown in Fig. 7 (dots) for j3 = 1/1.1. Here 
we define 

SO(t) = 0(t) - (£>>, (25) 

as the fluctuation from the mean. As the couplings are increased near their optimum 
values, the ringing disappears, and decorrelation times become better defined quantities, 
indicated by the solid curve with Ki = k 2 = 16. The corresponding HMC autocorrelation 
functions are shown as well. Parenthetically, this type of ringing can also occur in HMC 
simulations if one uses a constant trajectory length. The comparison to our HMC runs 
at this (3 is shown in Fig. 8. 

A comparison of auto-correlation functions for the total lattice magnetization M for 
global demons to the various implementations of HMC are shown in Fig. 9 for a selection 
of temperatures. It is clear that HMC does not significantly out perform global demons in 
terms of decorrelation, no matter how 'optimal' the trajectory length. The points in these 
curves indicate the actual number of data points. Hence while the number of global demon 
measurements is given by t, the optimal HMC runs have between one and two orders of 
magnitude smaller sampling rate in order to have similar decorrelation behavior. 



3.7 Critical Exponents 



The integrated autocorrelation times r are defined for a given quantity O as 



to —T 



1 " (5O(0)8O(t)) 

2 (50(0)50(0)} 



(26) 



Note that we are measuring r in units of total time evolved rather than number of trajec- 
tories. In Fig. 10 we present a ln-ln plot of the decorrelation times r vs. the correlation 
length £ for several observables. The fit parameters are indicated in Table 2, and are 
only indicative of the critical behavior, since they will depend strongly on systematic ef- 
fects. What is seen is that in almost every case, the global demon prefactor and critical 
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exponent are smaller than the HMC results. The integrated autocorrelation times are 
tabulated in Table 3. In Ref. it was argued that the critical exponent is unity when 
T is proportional to the lowest frequency mode in the system for a free field theory. The 
results in Tables 2-3 seem to be good evidence that their results are qualitatively correct 
in an interacting field theory as well. 

The critical behavior of global demon dynamics will also depend on the coupling 
strengths. In Fig. 11 we plot r as a function of k for simulations at several values of /3, 
on a 16 2 lattice. As k increases, the system tend to decorrelate faster, again generally 
saturating above k ~ ^/2n. In Fig. 12, the (3 dependence of simulations at k — 1, 4, 16, 64 
are indicated. The k = 1 runs have the highest decorrelation times as expected, but 
we also observe that the high temperature phase is rather insensitive to the value of the 
coupling. Although convergence of the trajectory to the correct ensemble will always 
depend strongly on the coupling strength, the decorrelation times of both weakly ergodic 
and strongly ergodic trajectories are very similar. The effects of critical slowing down are 
particularly noticeable in tq and te near j3 ~ 1. The peak in the te is closely related to 
the dip in AS in Fig. 4. The reason is that the diffusiveness measures the maximum rate 
at which the total energy S can change, so if AS <C 1, the potential energy E will change 
slowly. 

One might conclude from Figs. 11-12 that the coupling strength dependence is not too 
important, and that k ~ 2 is roughly equivalent to 128. Clearly, while the decorrelation 
times are indicative of the dynamics, they do not provide the complete picture of the 
situation. For example, information such as the ringing in the autocorrelation functions 
(see Fig. 7) average out, and are not strongly reflected in the value of r. We also see 
that simulations with very similar decorrelation times can have disparate values of the 
diffusiveness AS. But in all these guides, the tuning is consistently optimal for k ~ \/2n. 



4. Conclusions 



We have studied the global demon dynamical approach to simulating lattice regularized 
field theories. This method breaks away from the conventional Hamiltonian wisdom, 
defining a deterministically chaotic, time reversal invariant 'canonical' dynamics which 
rapidly fills configuration space with the desired ensemble. We have taken a particularly 
simple implementation of global demons using two coupling functions and examined its 
critical behavior, comparing to HMC using various trajectory lengths. 

We have found that the algorithm is very stable under tuning of the various parameters. 
In particular, once k is large enough, the quality of the results seem to be independent 
of k until the simulation becomes unstable. It appears that k ~ \^2n is a good rule of 
thumb for which the simulation will perform well in all regimes. In addition we found 
that the systematic errors were small (a few standard deviations), and that the critical 
slowing down properties of the algorithm were competitive with or better than the best 
implementations of HMC. The main fault with the algorithm is that it is not exact, i.e. 
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there is a systematic error associated with the numerical integration. This problem is 
addressed in the appendix. 

One advantage of this approach is that there is no barrier in principle to obtaining 
z < 1 (note our estimated critical exponents for S, E and Q in Table 2). In contrast, local 
algorithms such as HMC and HMD are limited by z ~ lH. In addition, the dynamical 
nature of the algorithm has allowed extensions to non-equilibrium situations [|ll|]. One 
possible improvement which we have not examined is Fourier acceleration. This technique 
improves z to for HMC in free field theory, and may help our approach. 

The numerical computations in this work were performed on the Cray Y-MP8 / 864 at 
the Ohio Supercomputer Center. We thank Aurel Bulgac, Robert Edwards, Rajan Gupta, 
Bill Hoover, Tony Kennedy, Greg Kilcup, Klaus Pinn, Junko Shigemitsu and Beth Thacker 
for useful conversations. This work was supported under DOE grants DE-AC02-ER01545 
and DE-FG02-91ER40608. 
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Appendix: Exact Dynamics? 



One major weakness of the global demon algorithm presented in this paper is the 
systematic error associated with finite step size. In HMC, this error is eliminated by per- 
forming a global metropolis hit before every measurement. Unfortunately, this technique 
has no hope of working for an algorithm with global demons, since the trajectories do 
not conserve energy. One might hope to use the pseudo-energy 8 in this capacity, since 
it is conserved. However, the memory term in Eq. (|13|) precludes this. For Hamiltonian 
systems, one can implement symplectic integrators PDA to render the dynamics exact, or 
one can introduce a global Metropolis hit as in HMC. For our non-Hamiltonian 'canonical' 
dynamics, this procedure is not so clear. Although we have been unable to satisfactorily 
solve this problem, we mention here one approach which we have tried. We present this 
method because, although it is presently numerically impractical, it appears to be correct 
in principle. In addition, it should be applicable to any dynamical simulation which does 
not conserve energy, and is quite different from the traditional method of using metropolis 
hits to ensure exactness. 

Consider the exact (i.e. dt = 0) equations of motion, <fi = F (<[>), which, by construction, 
preserve the measure p((f>) = exp(— (3f(4>)). When we discretize time, p is no longer 
preserved exactly. Let us assume that there exists some measure, p ^ p, which is preserved 
by the discretized equations of motion = M((j) n ), where M is the time evolution 
operator. We can now define a correction factor a such that, up to normalization, 

p(<P) = a(4>)p(<j>). (A.l) 

If we know a, we can obtain the exact expectation values of observables through convo- 
lution: 

!_Opd4 KOam = (Ooy 

[ ' fpd<f) Japd(f) (a)' ' 1 ' ' 

where the prime indicates evaluation in the p ensemble. The discretized Liouville (con- 
tinuity) equation tells us that after one integration time step, we preserve the measure 
P- 

Pn+1 d<p n+1 = p n d<p n . (A. 3) 

We can now use (A.l) and (A3) to solve for the correction factor a, through which we 
can measure exact observables: 

d(f) n+1 p n+1 fdM\ p n+1 

«n+l = "n— 77 = «n det -7— , (A.4) 

d<p n Pn \0<P J Pn 

where det(dM/d(f>) is just the Jacobian of the map <p n — > (j) n +i- Thus, by choosing a = 1 
at the initial point of the trajectory, we can compute a n at all subsequent points along 
the trajectory. 

Note that a should play a role very similar to the acceptance rate in HMC. When dt 
is small, the distribution p is very close to p, so a ~ 1 and all of the configurations will be 
of approximately equal statistical weight. This corresponds to high acceptance rates in 
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HMC When dt is large, on the other hand, a will be large in some regions and small in 
others, meaning that the simulation spends significant amounts of time in regions of low 
statistical importance. This is similar to the rejection of many proposed configurations 
when the acceptance is low. Also note that, although in general it is quite difficult, for our 
leapfrog algorithm the calculation of det(<9M/<90) is straight forward and only 0(volume). 
The correction factor can also be written as 



ot n = a a exp < in 

{ 8=0 



n-l 

/ ii i - , ■< i i t i i 

PHi + pHo) (A.5) 



dq i+1 dp i+1 



a 



dqi dpi 
exp{-/3 [£ n -£ + 0(At 2 )]}, 



where £ n is the pseudo-energy (|13|) evaluated at <p n . 

We have implemented this algorithm for small harmonic and anharmonic oscillator 
systems. In all runs, we found that, along the trajectory, the magnitude of a would 
occasionally rapidly decrease several orders of magnitude and never recover. We believe 
that this is due to an instability in the equations of motion: if the demons become too 
large at fixed k and dt, they begin to grow in an unbounded manner. It is quite likely 
that this means that a non-zero p does not exist. We have tried decreasing dt and various 
improvements in the equations of motion, but all of these fixes only decrease the frequency 
of drops - they do not eliminate them. 

In the scheme described above, it is not obvious how to perform a momentum refresh. 
This is because the absolute magnitude of a is not known, only its ratio to the previous 
alpha along the trajectory. We will now show (again, assuming the existence of p) that 
the correct procedure is to reset a to 1 after a gaussian momentum refresh. Assume, 
given p and a refreshing scheme, that the probability of a trajectory beginning at cf) is 
P{4>). In addition, let a(<fi) be the correctly normalized function a discussed above and 
a(0|0o) = «(0)/q;(0o) be the value of alpha obtained using (A 4) along the trajectory 
from cf) to 4> (note that a((p \(f) ) = 1). Then 

((D) = J dfoPjfo) f dtOpWaWfa) = N(Oa)> 

[ 1 J #o-P(0o) / d<f)p(<f))a((f)\<f)o) N(a)> 1 ' } 

where N — j d<f>o(P '(<po) j 'a(4>o)) . We attempted to test this scheme on the systems dis- 
cussed above, but we found that the statistical errors in the exact demon simulation 
tended to be larger than the systematic errors in a similar inexact simulation. Thus, it is 
presently more efficient to eliminate dt errors by extrapolation techniques than by using 
the exact algorithm. 



16 



References 

[1] S.Duane, A.D.Kennedy, B. J.Pendelton and D.Roweth, Phys. Lett. B195 (1987) 216. 
[2] P.B.Mackenzie, Phys. Lett. B226 (1989) 369. 

[3] M.Creutz and A.Gocksch, Phys. Rev. Lett. 63 (1989) 9; R.Gupta, G.Kilcup 
and S.Sharpe, Phys. Rev. D38 (1988) 1278; S.Gupta, A.Irback, F.Karsch and 
B.Petersson, Phys. Lett. B242 (1990) 437. 

[4] A.D.Kennedy and B. Pendelton, Nucl. Phys. (Proc. Suppl.) B20 (1991) 118. 

[5] G. Bathas and H. Neuberger, Phys. Rev. D 45 (1992) 3880. 

[6] J.Sloan, D.Kusnezov and A.Bulgac, (submitted) (1992). 

[7] A.Bulgac and D.Kusnezov, Phys. Rev. A42 (1990) 5045; D.Kusnezov, A.Bulgac and 
W.Bauer, Ann. Phys. (NY) 204 (1990) 155. 

[8] D.Kusnezov, Phys. Lett. 167A (1992). 

[9] D.Kusnezov, Phys. Lett. B (1992) in press; D.Kusnezov and A.Bulgac, Ann. Phys. 
(NY) 214 (1992) 180. 

[10] G.Parisi and Wu Yongshi, Sci. Sin. 24 (1981) 483; G.Parisi, Nucl. Phys. B180 (1981) 
378; Nucl. Phys. B205 (1982) 337. 

[11] W.G.Hoover, Computational Statistical Mechanics, (Elsevier, New York, 1991); 
S.Nose, Prog. Th. Phys. Suppl. 103 (1991) 1; A.Bulgac and D.Kusnezov, Phys. Lett. 
A151 (1990) 122. 

[12] The identification of stable periodic orbits is less clear. In practice such situations 
have not been encountered. In any case, this can be repaired by a periodic momentum 
refresh. 

[13] H.A.Posch and W.G.Hoover, Phys. Rev. A38, (1988) 473; Phys. Rev. A39, (1989) 
2175. 

[14] J.M.Kosterlitz and D.J.Thouless, J. Phys. C6 (1973) 1181; J.M.Kosterlitz, J. Phys. 
C7 (1974) 1046. 

[15] J. F. Fernandez, M. F. Ferreira and J. Stankiewicz, Phys. Rev. B34 (1986) 292; J. 
Kogut and J. Polonyi, Nucl. Phys. B265 (1986) 313; J. Tobochnik and G. V. Chester, 
Phys. Rev. B20 (1979) 3761; S.Miyashita, H.Nishimori, A.Kuroda and M. Suzuki, 
Prog, of Th. Phys. 60 (1978) 1669; S.Gupta, Nucl. Phys. B (1992). 

[16] R. Gupta, J. DeLapp and G. Batrouni, Phys. Rev. Lett. 61 (1988) 1996; W.Janke 
and K.Nather, Freie Universitat Berlin preprint FUB-HEP 7/91 (1991). 

[17] R. Gupta and C. Baillie, Los Alamos National Lab. preprint LA-UR-91-2364 (1991). 



17 



[18] The problem of ergodicity for weak and strong coupling limits has been addressed 
by K.Cho and J.D. Joannopoulos, Phys. Rev. A45 (1992) 7089, for Lennard- Jones 
potential systems. 

[19] S. Hands and R.J.Wensley, Phys. Rev. Lett. 63 (1989) 2169. 

[20] H. Yoshida, Phys. Lett. A150 (1990) 262. 



18 



TABLE 1. Comparison of observables between algorithms. For each f3, the four rows 
correspond to global demons, HMC-1, HMC-S, HMC-L, respectively. Measurements were 
performed on a 64 2 lattice, with 160K/T statistics. 



P E C v Q x 100 Xo x 100 | M | X m 



0.70 0.8287(1) 0.760(3) 6.241(2) 4.49(2) 0.0484(2) 12.2(1) 

0.8299(1) 0.765(5) 6.235(2) 4.50(3) 0.0487(3) 12.4(2) 

0.8297(2) 0.750(7) 6.233(2) 4.40(3) 0.0493(3) 12.7(1) 

0.8300(3) 0.686(14) 6.231(5) 4.36(7) 0.0490(3) 12.5(1) 



0.78 0.9573(1) 0.992(4) 4.725(2) 3.92(2) 0.0683(3) 24.3(2) 

0.9577(2) 0.994(8) 4.730(3) 3.93(3) 0.0683(6) 24.3(4) 

0.9573(3) 1.010(12) 4.731(3) 3.94(4) 0.0674(4) 23.6(3) 

0.9582(5) 1.061(30) 4.722(6) 4.04(9) 0.0683(5) 24.3(3) 



0.82 1.0243(1) 1.138(5) 3.981(2) 3.62(2) 0.0848(5) 37.2(4) 

1.0238(2) 1.13(1) 3.995(3) 3.61(3) 0.0828(9) 35.8(8) 

1.0244(3) 1.12(2) 3.987(3) 3.61(4) 0.0828(7) 35.7(6) 

1.0234(6) 1.09(4) 3.996(7) 3.53(10) 0.0837(7) 36.4(6) 



1/1.1 1.1757(2) 1.396(7) 2.431(2) 2.73(2) 0.162(1) 133(2) 

1.1750(3) 1.40(2) 2.448(4) 2.76(4) 0.151(3) 116(4) 

1.1756(5) 1.44(4) 2.443(5) 2.80(5) 0.161(2) 132(3) 

1.1777(10) 1.38(6) 2.418(10) 2.67(11) 0.166(2) 140(3) 



1/1.04 1.2632(3) 1.512(9) 1.638(3) 2.14(2) 0.295(3) 409(7) 

1.2618(5) 1.50(3) 1.659(5) 2.14(4) 0.276(7) 370(20) 

1.2618(7) 1.50(6) 1.657(6) 2.14(6) 0.281(4) 374(10) 

1.2630(13) 1.76(17) 1.647(12) 2.33(15) 0.291(5) 401(11) 
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TABLE 2. Comparison of estimated critical exponents z and prefactor A for global 
demons, HMC-1, HMC-S and HMC-L, where r = A£ z . Measurements were performed on 
a 64 2 lattice, with 160K/T statistics. 



A 



Lattice Magnetization M : 



Global Demons 1.3 2.4 

HMC-1 2.0 5.7 

HMC-S 1.3 6.9 

HMC-L 1.05 5.8 



Single Spin S : 

Global Demons 0.8 1.0 

HMC-1 1.5 1.6 

HMC-S 1.0 1.6 

HMC-L 1.0 4.3 



Potential Energy E : 

Global Demons 0.5 0.8 

HMC-1 1.0 1.9 

HMC-S 1.4 2.1 

HMC-L 1.3 14 



Topological Charge Q : 

Global Demons 0.8 0.8 

HMC-1 1.0 1.4 

HMC-S 1.3 1.6 

HMC-L 1.3 8.2 
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TABLE 3. Auto-correlation times for total magnetization M, a single spin S, topological 
charge Q and internal energy E, and the magnetization correlation length £. Measure- 
ments were performed on a 64 2 lattice, with 160K/T statistics. 
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FIGURES 



1. Time evolution of 10 3 connected points evolving according to 'canonical' rather than 
microcanonical dynamics for (a) the 1-d harmonic oscillator H = (p 2 + <? 2 )/2 (left 
column; t — points lie on the unit circle), and (b) the 577(2) Hamiltonian H = J 2 /2 
(right column; t — points lie on the circle J z = 0.5), both at f3 — 1. The phase 
space in (a) is (q,p) and in (b) it is the sphere parameterized by (J x , J y , J z ). The 
rapid spreading of neighboring points is characteristic of global demon dynamics. 
The characteristic time scale is t ~ 1. 

2. Finite time step extrapolations, demonstrating the global 0(dt 2 ) leapfrog error for 
(a) potential energy, (b) topological charge and (c) magnetization. 

3. Diffusiveness of a global demon trajectory versus the coupling strength, for selected 
values of (3. AS is measured at intervals of T = 1. The 'optimal' value of \J (AS 2 ) = 
a is indicated by the dashed line. As can be seen, the dynamics is not strongly 
affected by critical slowing down. Saturation occurs when n a ~ 0(y/2n) ~ 23. The 
optimal coupling can be seen to be independent of both (3 and the phase transition. 

4. Diffusiveness of the global demons versus f3. The dip at (3 ~ 1 is a result of the 
Kosterlitz-Thouless phase transition and critical slowing down. The dynamics is 
not strongly effected by the transition, so no additional tuning is required, and k 
can be taken as fixed for all (3. 

5. Measurement dependence on the coupling k at (3 = 1/1.1. 

6. Diffusiveness of a global demon trajectory a function of the trajectory length T, 
at (3 — 1, Ki — k 2 — 16. Saturation can be seen, indicating an optimal trajectory 
length of T ~ 0(1). 

7. Autocorrelation functions for potential energy, topological charge, magnetization 
and spin, at (3 = 1/1.1 with k± — k 2 — 1 (dots) and K\ — k 2 — 16 (solid). As can 
be seen, the ringing vanishes as the coupling increases. 

8. Autocorrelation functions for potential energy, topological charge, magnetization 
and spin, at (3 = 1/1.1 for global demons with K\ — k 2 — 16 (solid), HMC-1 
(dots), HMC-S (crosses) and HMC-L (boxes). For global demons and HMC-1, 
measurements are made every t — 1, while for HMC-S, L, the boxes and crosses 
indicate the actual number of data points. 

9. Lattice magnetization auto-correlation function at selected temperatures for global 
demons (solid), HMC-1 (dots), HMC-S (dashes) and HMC-L (boxes). The time axis 
for (3 = 0.7 has been scaled by a factor of 0.1 to magnify the short time behavior. 

10. The behavior of r = a£ z is plotted near the phase transition on a 64 2 lattice for 
(a) total magnetization, (b) the spin at the origin, (c) topological charge and (d) 
potential energy. In each figure, we indicate the results for HMC (boxes), HMC— s 
(diamonds), HMC—/ (crosses) and global demons (circles). Critical exponents can 
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be extracted from a linear fit. The results of the numerical fits are given in Table 
2. Results are on a 64 2 lattice with 1Q0K/T statistics. 



11. Decorrelation times versus coupling strength for (3 = 0.5 (squares), 0.7 (diagonal 
crosses), 1/1.1 (diamonds), 2.0 (vertical crosses) and 4.0 (circles) for (a) total mag- 
netization, (b) a single spin, (c) topological charge and (d) potential energy. 

12. Decorrelation times versus (3 for couplings k± — k 2 — 1 (crosses), 4 (diamonds), 16 
(squares) and 64 (circles) for (a) total magnetization, (b) a single spin, (c) topolog- 
ical charge and (d) potential energy. The rise at /3 ~ 1 is critical slowing down, and 
is especially evident for the energy and topological charge. 
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